home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
SCIENTIF
/
H381.ZIP
/
GSRC208A.ZIP
/
GAUSSW.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-26
|
19KB
|
600 lines
#include "copyleft.h"
/*
GEPASI - a simulator of metabolic pathways and other dynamical systems
Copyright (C) 1989, 1992 Pedro Mendes
*/
/*************************************/
/* */
/* GWSIM - Simulation */
/* MS-WINDOWS front end */
/* */
/* Gauss reduction */
/* dialog box */
/* */
/* QuickC/WIN 1.0 */
/* */
/* (include here compilers that */
/* compiled GWSIM successfully) */
/* */
/*************************************/
#include <windows.h>
#include <math.h>
#include <string.h>
#include "globals.h"
#include "gep2.h"
#include "strtbl.h"
#define ALMOST_ZERO 1.0e-7
GLOBALHANDLE hRSto; /* handle to memory block w/ rstoi */
GLOBALHANDLE hSto; /* handle to memory block w/ stoi */
GLOBALHANDLE hMl; /* handle to memory block w/ ml */
GLOBALHANDLE hLm; /* handle to memory block w/ lm */
GLOBALHANDLE hLD; /* handle to memory block w/ ld */
GLOBALHANDLE hMetn; /* handle to memory block w/ metn */
GLOBALHANDLE hStepn; /* handle to memory block w/ stepn */
GLOBALHANDLE hStepst; /* handle to memory block w/ stepst */
GLOBALHANDLE hImet; /* handle to memory block w/ imet */
GLOBALHANDLE hRs; /* handle to memory block w/ rs */
GLOBALHANDLE hLoo; /* handle to memory block w/ loo */
float huge *rstoi; /* pointer to work with rstoi array */
float huge *ml; /* pointer to work with ml array */
float huge *lm; /* pointer to work with lm array */
float huge *ld; /* pointer to work with ld array */
int huge *imet; /* pointer to work with imet vector */
int huge *loo; /* pointer to work with loo array */
int (huge *rs)[MAX_STEP][MAX_MOL]; /* pointer to work with rs array */
char (huge *metn)[NAME_L]; /* pointer to work with metname array */
char (huge *stepn)[NAME_L]; /* pointer to work with metname array */
char (huge *stepst)[256]; /* pointer to work with stepst array */
int ur[MAX_MET]; /* permut. on metabolites u -> g */
int uc[MAX_MET]; /* permutations on steps u -> g */
void rowsw( int r1, int r2, int c );
void colsw( int c1, int c2, int r );
void rowsw_m( int r1, int r2, int c );
void colsw_m( int c1, int c2, int r );
int emptyrow( int rw, int nsteps);
void invml11( void );
void calc_ld( void );
void init_moiety( void );
void lindep( void );
int gauss( void );
void more_ext( void);
void int_ord( void );
void ext_ord( void );
void initreds( void );
void virt_step( void );
int reduce( int swapnames );
#pragma alloc_text( CODE14, rowsw, colsw, rowsw_m, colsw_m, emptyrow, invml11, calc_ld, init_moiety, lindep, gauss, int_ord, ext_ord, more_ext, initreds, virt_step, reduce )
void rowsw( int r1, int r2, int c )
{
register int j;
float dummy;
double dumb;
int dum;
for(j=0;j<c;j++)
{
dum = stoi[r1*MAX_MET + j]; /* switch rows in stoi and */
stoi[r1*MAX_MET + j] = stoi[r2*MAX_MET + j];
stoi[r2*MAX_MET + j] = dum;
dummy = rstoi[r1*MAX_MET + j]; /* switch rows in rstoi & */
rstoi[r1*MAX_MET + j] = rstoi[r2*MAX_MET + j];
rstoi[r2*MAX_MET + j] = dummy;
}
j = ur[r1]; /* switch rows in ur and */
ur[r1] = ur[r2];
ur[r2] = j;
}
void colsw( int c1, int c2, int r )
{
register int j;
float dummy;
int dum;
for( j=0; j<r; j++ )
{
dum = stoi[j*MAX_MET + c1]; /* switch cols in stoi and */
stoi[j*MAX_MET + c1] = stoi[j*MAX_MET + c2];
stoi[j*MAX_MET + c2] = dum;
dummy = rstoi[j*MAX_MET + c1]; /* switch cols in rstoi & */
rstoi[j*MAX_MET + c1] = rstoi[j*MAX_MET + c2];
rstoi[j*MAX_MET + c2] = dummy;
}
j = uc[c1]; /* switch cols in uc too */
uc[c1] = uc[c2];
uc[c2] = j;
}
void rowsw_m( int r1, int r2, int c ) /* switch rows in ml */
{
register int j;
float dummy;
for(j=0;j<c;j++)
{
dummy = ml[r1*MAX_MET + j];
ml[r1*MAX_MET + j] = ml[r2*MAX_MET + j];
ml[r2*MAX_MET + j] = dummy;
}
}
void colsw_m( int c1, int c2, int r ) /* switch cols in ml */
{
register int j;
float dummy;
for(j=0;j<r;j++)
{
dummy = ml[j*MAX_MET + c1];
ml[j*MAX_MET + c1] = ml[j*MAX_MET + c2];
ml[j*MAX_MET + c2] = dummy;
}
}
int emptyrow( int rw, int nsteps)
{
register int j, ct; /* ct counts entries != 0 */
for( ct=0, j=0; j<nsteps; j++)
if ( fabs( rstoi[rw*MAX_MET + j] ) > ALMOST_ZERO ) ct++;
return ( ct );
}
void invml11( void )
{ /* as ml is lower triangul.*/
register int i,j,k; /* inverting is simply to */
float acum; /* forward-substitute with */
/* the identity matrix */
for(i=0;i<nmetab; i++) ml[i*MAX_MET + i] = 1 ;
for( j=0; j<indmet; j++)
for( k=0; k<indmet; k++)
{
for( acum=0, i=0; i<k; i++ )
acum += ml[k*MAX_MET + i] * lm[i*MAX_MET + j];
lm[k*MAX_MET + j] = ( (k==j) ? (float) 1.0 : - acum ) / ml[k*MAX_MET + k];
}
}
void calc_ld( void )
{
register int i,j,k;
for( i=indmet; i<nmetab; i++ ) /* first calculate L0 */
for( j=0; j<indmet; j++ ) /* which is */
{
ld[i*MAX_MET + j] = 0;
for( k=0; k<indmet; k++ ) /* -1 */
ld[i*MAX_MET + j] += ml[i*MAX_MET + k] * lm[k*MAX_MET + j]; /* L0 = L21 * (L11) */
}
for( i=indmet; i<nmetab; i++ ) /* Now map the lin. dep. */
{
for( j=0; j<indmet; j++ ) /* which is */
ld[i*MAX_MET + j] = -ld[i*MAX_MET + j]; /* -L0 */
for( j=indmet; j<nmetab; j++ ) /* concateneted with the */
ld[i*MAX_MET + j] = (i == j)
? (float) 1.0 : (float) 0.0; /* m0 -> m part of Id. */
}
}
void init_moiety( void )
{
register int i,j;
for( i=indmet; i<nmetab; i++)
{
moiety[i] = (float) 0.0;
for( j=0; j<nmetab; j++ )
moiety[i] += ld[i*MAX_MET + j] * xu[j];
}
}
void lindep( void )
{
invml11(); /* invert ml11 into lm */
calc_ld(); /* calculate L2 * -L1' */
init_moiety(); /* setup the moiety conc. */
}
int gauss( void )
{
register int i, j;
int k, flag;
float m;
for( k=0; k<nsteps+1; k++)
{
for( flag=0, i=k; i<nmetab; i++ )
if( fabs(rstoi[i*MAX_MET + k]) > ALMOST_ZERO )
{
flag = 1;
break; /* suitable divisor */
}
if ( flag )
{
if ( i != k )
{
rowsw( k, i, nsteps ); /* switch row k with i */
rowsw_m( k, i, nmetab );
}
}
else
{
do
{
for( j=k+1; j<nsteps; j++ )
if( fabs(rstoi[k*MAX_MET + j]) > ALMOST_ZERO )
{
flag = 1;
break; /* suitable value */
}
if ( flag && ( j != k ) )
{
colsw( j, k, nmetab ); /* switch col j with k */
colsw_m( j, k, nmetab );
}
if( !flag )
{
for( i=0; i<nmetab; i++ ) /* get rid of empty rows */
{
if ( ! emptyrow( i, nsteps ) ) /* all entries in row = 0 */
{
for( j=i+1; j<nmetab; j++ )
if ( emptyrow( j, nsteps ) ) /* j = first non-empty row */
{
rowsw( i, j, nsteps ); /* switch row i with row j */
rowsw_m( i, j, nmetab );
flag = 2;
break;
}
}
}
}
if ( !flag )
{
indmet = k; /* # of independent metabs */
return(-1); /* flag conservation & ret */
}
}
while( flag != 1 );
}
for( i=k+1; i<nmetab; i++)
{
if ( fabs(rstoi[k*MAX_MET + k]) > ALMOST_ZERO )
{
m = rstoi[i*MAX_MET + k] / rstoi[k*MAX_MET + k]; /* calculate the divisor */
ml[i*MAX_MET + k] = m; /* and keep its symmetric */
if( m != (float) 0.0 )
{
for( j=k; j<nsteps; j++ )
{
rstoi[i*MAX_MET + j] = rstoi[i*MAX_MET + j] - m*rstoi[k*MAX_MET + j]; /* reduce another row */
if( fabs(rstoi[i*MAX_MET + j]) <= ALMOST_ZERO ) /* then make it a true zero*/
rstoi[i*MAX_MET + j] = (float) 0.0;
}
}
}
}
}
return 0; /* job done: no lin. dep.! */
}
void int_ord( void ) /* start by putting all internal */
{ /* metabolites as the first */
register int i,j,k; /* m rows in stoi */
for( i=0, k=0; i<totmet; i++ )
if (intmet[i])
{
ur[k] = i; /* record the old number */
for( j=0; j<nsteps; j++ )
stoi[k*MAX_MET + j] = stoiu[i*MAX_MET + j];
k++;
}
nmetab = k; /* the no.of internal metabolites */
}
void ext_ord( void ) /* now put the declared exernal */
{ /* metabolites in the last rows */
register int i,j,k; /* of stoi */
k = nmetab;
for( i=0; i<totmet; i++ )
if (!intmet[i])
{
ur[k] = i; /* record the old number */
for( j=0; j<nsteps; j++ )
stoi[k*MAX_MET + j] = stoiu[i*MAX_MET + j];
k++;
}
nextmet = totmet - nmetab; /* the no.of external metabolites */
}
void initreds( void )
{
register int i,j;
for( i=0; i<nmetab; i++)
{
for( j=0; j<nsteps; j++)
rstoi[i*MAX_MET + j] = (float) stoi[i*MAX_MET + j]; /* init rstoi from stoi */
for( j=0; j<nmetab; j++)
ml[i*MAX_MET + j] = (float) 0; /* init matrix of multipl. */
}
}
void more_ext( void)
{
register int i;
for( i=0; i<nmetab; i++ ) /* get rid of empty rows */
{ /* which are xtrnl metabs. */
if ( ! emptyrow( i, nsteps ) ) /* all entries in row = 0 */
{
rowsw( i, nmetab-1, nsteps ); /* put row at the end */
intmet[ur[nmetab]] = 0; /* signal this 1 is xtrnl! */
nmetab--; /* signal 1 less int. met. */
nextmet++; /* and 1 more extern. met. */
}
}
}
void virt_step( void )
{
register int i, j, ct;
for( i=0; i<nsteps; i++ ) /* get rid of empty cols */
{ /* which are 0 rate steps */
for( ct=0, j=0; j<nmetab; j++)
if ( fabs( rstoi[j*MAX_MET + i] ) > ALMOST_ZERO ) ct++; /* ct counts entries != 0 */
if ( !ct ) /* all entries in col = 0 */
{
colsw( i, nsteps-1, nmetab ); /* put col at the end */
nsteps--; /* one less active step */
}
}
}
int reduce( int swapnames )
{
int i, j, k, ix, s;
int kinetype[MAX_STEP]; /* type of kinetics (user numb.) */
/* first lets allocate memory for mirrors and aux matrices */
hRSto = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * MAX_STEP * sizeof( float ) );
if( hRSto == NULL )
return IDS_ERR_OUT_OF_MEM;
hMetn = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * NAME_L * sizeof( char ) );
if( hMetn == NULL )
{
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hStepn = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_STEP * NAME_L * sizeof( char ) );
if( hStepn == NULL )
{
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hStepst = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_STEP * 256 * sizeof( char ) );
if( hStepst == NULL )
{
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hMl = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * MAX_MET * sizeof( float ) );
if( hMl == NULL )
{
GlobalFree( hStepn );
GlobalFree( hStepst );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hLm = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * MAX_MET * sizeof( float ) );
if( hLm == NULL )
{
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hLD = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * MAX_MET * sizeof( float ) );
if( hLD == NULL )
{
GlobalFree( hLm );
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hImet = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_MET * sizeof( int ) );
if( hImet == NULL )
{
GlobalFree( hLD );
GlobalFree( hLm );
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hRs = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_STEP * MAX_MOL * sizeof( int ) );
if( hRs == NULL )
{
GlobalFree( hImet );
GlobalFree( hLD );
GlobalFree( hLm );
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
hLoo = GlobalAlloc( GMEM_MOVEABLE | GMEM_ZEROINIT, (DWORD) MAX_STEP * MAX_MET * sizeof( int ) );
if( hLoo == NULL )
{
GlobalFree( hRs );
GlobalFree( hImet );
GlobalFree( hLD );
GlobalFree( hLm );
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return IDS_ERR_OUT_OF_MEM;
}
/* and lock the memory blocks */
rstoi = (float huge *) GlobalLock( hRSto );
metn = (char (huge *)[NAME_L]) GlobalLock( hMetn );
stepn = (char (huge *)[NAME_L]) GlobalLock( hStepn );
stepst = (char (huge *)[256]) GlobalLock( hStepst );
ml = (float huge *) GlobalLock( hMl );
lm = (float huge *) GlobalLock( hLm );
ld = (float huge *) GlobalLock( hLD );
imet = (int huge *) GlobalLock( hImet );
rs = (int (huge *)[MAX_MET][MAX_MOL]) GlobalLock( hRs );
loo = (int huge *) GlobalLock( hLoo );
for( i=0; i<nsteps; i++) uc[i] = i; /* setup reaction perm. vector */
for( i=0; i<totmet; i++) ur[i] = i; /* setup metabolite perm. vector */
int_ord(); /* put int. metabs. at the beggin.*/
ext_ord(); /* put ext. metabs. at the end */
initreds(); /* reset rstoi and ml */
more_ext(); /* spot implicit external metabs. */
virt_step(); /* spot virtual steps */
gauss(); /* stoi -> rstoi by gaussian red. */
/* make the mirrors be indexed by the new mapping */
/* if swapnames remap metabolite and step names */
for( i=0; i<nsteps; i++ )
{
if( swapnames )
{
lstrcpy( (LPSTR) stepn[i], (LPSTR) stepname[uc[i]] );
lstrcpy( (LPSTR) stepst[i], (LPSTR) stepstr[uc[i]] );
kinetype[i] = kinetu[uc[i]];
for(j=0;j<totmet;j++)
loo[MAX_STEP*i+j] = (*loop)[uc[i]][ur[j]];
for(j=0;j<MAX_MOL;j++)
{
ix = (*rstr)[uc[i]][j];
if( ix == 0 )
{
(*rs)[i][j] = 0;
continue;
}
s = ix/ abs( ix ); /*the signal*/
ix = abs( ix )-1; /*the index */
for( k=0; k<totmet; k++ )
if( ix == ur[k] ) break;
(*rs)[i][j] = s*(k+1);
}
}
else
{
lstrcpy( (LPSTR) stepn[i], (LPSTR) stepname[i] );
lstrcpy( (LPSTR) stepst[i], (LPSTR) stepstr[i] );
kinetype[i] = kinetu[i];
for(j=0;j<totmet;j++)
loo[MAX_STEP*i+j] = (*loop)[i][j];
for(j=0;j<MAX_MOL;j++)
(*rs)[i][j] = (*rstr)[i][j];
}
}
/* if( swapnames )
{
for( i=0; i<nsteps; i++)
for(j=0;j<MAX_MOL;j++)
{
if( (*rstr)[uc[i]][j] > 0 )
(*rs)[i][j] = ur[(*rstr)[uc[i]][j] - 1] + 1;
else
if( (*rstr)[uc[i]][j] < 0 )
(*rs)[i][j] = -ur[-(*rstr)[uc[i]][j] - 1] - 1;
else
(*rs)[i][j] = 0;
}
} */
for( i=0; i<totmet; i++ )
if( swapnames )
{
lstrcpy( (LPSTR) metn[i], (LPSTR) metname[ur[i]] );
imet[i] = intmet[ur[i]];
}
else
{
lstrcpy( (LPSTR) metn[i], (LPSTR) metname[i] );
imet[i] = intmet[i];
}
/* copy everything back to the original structures */
_fmemcpy( (void __far *) metname,
(void __far *) metn,
(size_t) MAX_MET * NAME_L * sizeof( char ) );
_fmemcpy( (void __far *) stepname,
(void __far *) stepn,
(size_t) MAX_STEP * NAME_L * sizeof( char ) );
_fmemcpy( (void __far *) stepstr,
(void __far *) stepst,
(size_t) MAX_STEP * 256 * sizeof( char ) );
_fmemcpy( (void __far *) kinetu,
(void __far *) kinetype,
(size_t) MAX_STEP * sizeof( int ) );
_fmemcpy( (void __far *) intmet,
(void __far *) imet,
(size_t) MAX_MET * sizeof( int ) );
_fmemcpy( (void __far *) rstr,
(void __far *) rs,
(size_t) MAX_STEP * MAX_MOL * sizeof( int ) );
for(i=0;i<nsteps;i++)
for(j=0;j<totmet;j++)
(*loop)[i][j] = (unsigned char) loo[MAX_STEP*i+j];
lindep(); /* work out linear dependencies */
/* unlock and free the mirrors */
GlobalUnlock( hLoo );
GlobalUnlock( hRs );
GlobalUnlock( hImet );
GlobalUnlock( hLD );
GlobalUnlock( hLm );
GlobalUnlock( hMl );
GlobalUnlock( hStepst );
GlobalUnlock( hStepn );
GlobalUnlock( hMetn );
GlobalUnlock( hRSto );
GlobalFree( hLoo );
GlobalFree( hRs );
GlobalFree( hImet );
GlobalFree( hLD );
GlobalFree( hLm );
GlobalFree( hMl );
GlobalFree( hStepst );
GlobalFree( hStepn );
GlobalFree( hMetn );
GlobalFree( hRSto );
return 0;
}